#define _targ_ess 100
template<class Approx,class APriori,class Likelihoodms>
class Probit_MSadap : public Density::GeomBridge
{
		friend class Density::GeomBridge;
	public:
		Probit_MSadap(double *Y, mat X, Distribution::Distribution *S, APriori *P, Likelihoodms *Lik,int ds) : Density::GeomBridge(Y, X, S){ 
			mat s=(*P).Get_s();
			_s=P->Get_s();
			_p=X.n_cols-1;
			_ds=ds;
			_d=new boost::math::normal_distribution<>(0,1);
			Set_p(X.n_cols-1);
			_is=inv(_s);
			_Lik=Lik;
			fL=new Approx;
			mat m(_p+1,1);
			ds=logfact(_ds);
			d=logfact(_p);
			m.fill(0);
			(*fL)(X,Y,m,_s(0,0),m);			
			_init=fL->Get_Mu();
			_Px=P;
		}
		~Probit_MSadap(){
			delete fL;
		}
		double Likelihood(mat gamma)
		{
			int p=Get_p();
			mat X=Get_X();
			mat gammat(p+1,1);
			mat h=sum(gamma,0);
			double res=0;
			int nb=h(0,0)+1;
			int n=X.n_rows;	
			mat X2(n,nb);
			mat init(nb,1);
//			cout << "test"<< p+1;
			int kk=1;
			for(int i=0;i<(p+1);i++)
			{
					
				if(i!=0){
					gammat(i,0)=gamma(i-1,0);
					if(gamma(i-1,0)!=0)
					{
						X2.col(kk)=X.col(i);
						init(kk,0)=_init(i,0);
						kk++;
					}

				}else{
					gammat(i,0)=1;
					X2.col(0)=X.col(0);
					init(0,0)=_init(0,0);
				}
			}
			//cout << gammat.t();
	
			mat Xt=X*diagmat(gammat);
			mat m(nb,1);
			mat s(nb,nb);
			s.fill(_s(0,0));
			s=diagmat(s);
			m.zeros();
			double *Y=Get_Y();

			(*fL)(X2,Y,m,_s(0,0),m);			
			mat Sig=(*fL).Get_Sig();	
			mat Mu=(*fL).Get_Mu();
		//	cout << Mu.t();

	
			//Distribution::Gaussian G(Xt.n_cols,m,s);
			Distribution::Gaussian G(nb,Mu,Sig);
			APriori F(nb,m,s);
			
			//cout << "test" << nb-prop.n_rows;
			double sum2=0;
			
		//	prop=(G.scrambled(_m)).t();
	//		cout << prop;

			mat prop=G.r(1);
			mat beta=add_mat(prop,gammat,p+1);
			sum2=Lik_Prob(beta,Xt)+F.d(prop,1)-G.d(prop,1);
			//	gammat.fill(1);
			vector<double> omega;
			omega.push_back(pow(exp(sum2),2));	
			double ess=0;
			int k=1;
			while(ess<_targ_ess)
			{
				prop=G.r(1);
				beta=add_mat(prop,gammat,p+1);
				double o=Lik_Prob(beta,Xt)+F.d(prop,1)-G.d(prop,1);
				double t=log_add(o,sum2);
				sum2=t;
				omega.push_back(pow(exp(o),2));
				mat normo(k,1);
				for(int j=0;j<k;j++)
				{
					normo(j,0)=omega[j]/exp(sum2);
				}
				ess=as_scalar(1/sum(normo,0));
				cout << " "<< ess;
				k++;
			}
			res=sum2-log(k);
			//cout << res << "\n";
			

			return res;

		}
		inline double Lik_Prob(mat theta1,mat Xt)
		{
			double* Y=Get_Y();
			double u=(*_Lik)(theta1,Xt,Y);
			/*double L=0;
			double sum=0;
			int m=theta1.n_cols;
			//Xb
			int n=Get_n();
			double* Y=Get_Y();
			mat foo(1,Xt.n_cols);
			for(int i=0;i<n;i++)
			{	
				foo= Xt.row(i);
				sum=dot(foo,theta1);
				if(Y[i]==1){
					L+=log(Phi(sum));
				}else{
					L+=log(1-Phi(sum));
				}	
			}
			//cout << exp(L);
			return L;
			*/
			return u;

		}
		inline double Phi(double x)
		{
			double t=cdf((*_d),x);
		//	double t=0.5*(1+Erfm(x/sqrt(x)));
			if(t==1)
			{
				t=0.999999;
			}
			if(t==0)
			{
				t=0.000001;
			}
			return t;
		}

		double Prior(mat theta){
		    
			double t=0;
			double m=0.5;
			double i=sum(sum(theta));
			if(i<=_ds)
			{
				double t=logfact(i)+(logfact(_ds-i))-ds+i*log(m)+(_ds-i)*log(1-m)+logfact(i)+(logfact(_p-i))-d;
				return t;
 			}else{
				return -numeric_limits<double>::infinity();

			}		
			
	//		return -0.693;
		}
	
		mat  GradLik(mat theta)
		{
			cout << "Grad ProbitMS deprecated";
			return theta;
		}

	private:
		Approx *fL;
		Likelihoodms *_Lik;
		boost::math::normal_distribution<> *_d;
		Distribution::Distribution *_Px;
		mat _s;
		int _ds;
		int ds;
		int d;
		mat _init;
		int _m;
		int _p;
		mat _is;
	
};
